APS/123-QED 



Crystal nucleation of hard spheres using molecular dynamics, umbrella sampling and 
forward flux sampling: A comparison of simulation techniques 

L. Filion, M. Hermes, R. Ni, and M. Dijkstra 
Soft Condensed Matter, Debye Institute for NanoMaterials Science, 
Utrecht University, Princetonplein 5, 3584 CC Utrecht, The Netherlands 
(Dated: June 16, 2010) 

Over the last number of years several simulation methods have been introduced to study rare 
events such as nucleation. In this paper we examine the crystal nucleation rate of hard spheres 
using three such numerical techniques: molecular dynamics, forward flux sampling and a Bennett- 
Chandler type theory where the nucleation barrier is determined using umbrella sampling simula- 
tions. The resulting nucleation rates are compared with the experimental rates of Harland and Van 
Megen [J. L. Harland and W. van Megen, Phys. Rev. E 55, 3054 (1997)], Sinn et al. [C. Sinn 
et al., Prog. Colloid Polym. Sci. 118, 266 (2001)] and Schatzel and Ackerson [K. Schatzel and 
B.J. Ackerson, Phys. Rev. E, 48, 3766 (1993)] and the predicted rates for monodisperse and 5% 
polydisperse hard spheres of Auer and Frenkel [S. Auer and D. Frenkel, Nature 409, 1020 (2001)]. 
When the rates are examined in long-time diffusion units, we flnd agreement between all the theo- 
retically predicted nucleation rates, however, the experimental results display a markedly different 
behaviour for low supersaturation. Additionally, we examined the pre-critical nuclei arising in the 
molecular dynamics, forward flux sampling, and umbrella sampling simulations. The structure of 
the nuclei appear independent of the simulation method, and in all cases, the nuclei contain on av- 
erage significantly more face-centered-cubic ordered particles than hexagonal-close-packed ordered 
particles. 



I. INTRODUCTION 

Nucleation processes are ubiquitous in both natural 
and artificially-synthesized systems. However, the occur- 
rence of a nucleation event is often rare and difficult to 
examine both experimentally and theoretically. 

Colloidal systems are almost ideal model systems for 
studying nucleation phenomena. Nucleation and the pro- 
ceeding crystallization in such systems often take place 
on experimentally accessible time scales, and due to the 
size of the particles, they are accessible to a wide variety 
of scattering and imaging techniques, such as (confocal) 
microscopy,'^ holography,El and light and x-ray scattering. 
Additionally, progress in particle synthesis,'' solvent ma- 
nipulation, and the application of external fields* allows 
for significant control over the interparticle interactions, 
allowing for the study of a large variety of nucleation 
processes. 

One such colloidal system is the experimental realiza- 
tion of "hard" spheres comprised of sterically stabilized 
polymethylmethacrylate (PMMA) particles suspended in 
a liquid mixture of decaline and carbon disulfide.'^' Exper- 
imentally, the phase behaviour of such a system has been 
examined by Pusey and Van MegerP and maps well onto 
the phase behaviour predicted for hard spheres. Specifi- 
cally when the effective volume fraction of their system is 
scaled to reproduce the freezing volume fraction of hard 
spheres (rj = 0.495) the resulting melting volume fraction 
is 77 = 0.545 ± 0.0032- which is in good agreement with 
that predicted for hard spheres.!^ The nucleation rates 
have been measured using light scattering by Harland 
and Van Megen,'^' Sinn et a/.P Schatzel and AckersorP 
and predicted theoretically by Auer and Frenkel.'^ 

On the theoretical side, hard-sphere systems are one of 



the simplest systems which can be applied to the study 
of colloidal and nanoparticle systems, and generally, to- 
wards the nucleation process itself. As such, it is an ideal 
system to examine various computational methods for 
studying nucleation, and comparing the results with ex- 
perimental data. Such methods include, but are not lim- 
ited to, molecular dynamics (MD) simulations, umbrella 
samphng (US), forward flux sampling (FFS), and tran- 
sition path sampling (TPS). It is worth noting here that 
Auer and Frenkel'^ used umbrella sampling simulations 
to study crystal nucleation of hard spheres and found a 
significant difference between their predicted rates and 
the experimental rates of Refs. HHSl However, it was un- 
clear where this difference originated. In this paper we 
compare the nucleation rates for the hard-sphere system 
from MD, US and FFS simulations with the experimen- 
tal results of Refs. HHSl We demonstrate that the three 
simulation techniques are consistent in their prediction 
of the nucleation rates, dispite the fact that they treat 
the dynamics differently. Thus we conclude that the dif- 
ference between the experimental and theoretical nucle- 
ation rates identified by Auer and Frenkel is not due to 
the simulation method. 

A nucleation event occurs when a statistical fluctua- 
tion in a supersaturated liquid results in the formation of 
a crystal nucleus large enough to grow out and continue 
crystallizing the surrounding fluid. In general, small crys- 
tal nuclei are continuously being formed and melting back 
in a liquid. However, while most of these small nuclei 
will quickly melt, in a supersaturated liquid a fraction of 
these nuclei will grow out. Classical nucleation theory is 
the simplest theory available for describing this process. 
In CNT it is assumed that the free energy for making a 
small nuclei is given by a surface free energy cost which 
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is proportional to the surface area of the nucleus and a 
bulk free energy gain proportional to its volume. More 
specifically, according to CNT the Gibbs free energy dif- 
ference between a homogeneous bulk fluid and a system 
containing a spherical nucleus of radius R is given by 

AG{R)^i7rjR^~^7T\An\p,R^ (1) 

where | A/i| is the difference in chemical potential between 
the fluid and solid phases, ps is the density of the solid, 
and 7 is the surface tension of the fluid-solid interface. 
This free energy difference is usually referred to as the 
nucleation barrier. From this expression, the radius of 
the critical cluster is found to be R* =27/ |A^| ps and 
the barrier height is AG* = l6TT-f^ /3p'j^ l^fJ-l^ ■ Note that 
there is no system si ze dep endence in CNT. 

Umbrella samplin^^iE^ is a method to examine the 
nucleation process from which the nucleation barrier is 
easily obtained. The predicted barrier can then be used 
in combination with kinetic Monte Carlo (KMC) or MD 
simulations to determine the nucleation rate.l^ In US an 
order parameter for the system is chosen and configura- 
tion averages for sequential values of the order parameter 
are taken. In order to facilitate such averaging, the sys- 
tem is biased towards particular regions in configuration 
space. The success of the method is expected to depend 
largely on the choice of order parameter and biasing po- 
tential. Note that the free energy barrier is only defined 
in equilibrium, and thus is only applicable to systems 
which are in (quasi-) equilibrium. 

Forward fiux samplin^l^Hls] jg rnethod of studying 
rare events, such as nucleation, in both equilibrium and 
non-equilibrium systems. Using FFS, the transition rate 
constants (eg. the nucleation rate) for rare events can be 
determined when brute force simulations are difficult or 
even not possible. In FFS, a reaction coordinate Q (sim- 
ilar to the order parameter in US) is introduced which 
follows the rare event. The transition rate between phase 
A and B is then expressed as a product of the fiux (<&aAo ) 
of trajectories crossing the A state boundary, typically 
denoted Xq, and the probability {P{Xb\Xq)) that a tra- 
jectory which has crossed this boundary will reach state 
B before returning to state A. Thus the transition rate 
constant is written as 

kAB^^AXoP{^B\Xo). (2) 

Forward flux sampling facilitates the calculation of prob- 
ability P{Xb\Xo) by breaking it up into a set of probabili- 
ties between sequential values of the reaction coordinate. 
Little information regarding the details of the nucleation 
process is required in advance, and the choice of reac- 
tion coordinate is expected to be less important than the 
order parameter in US. Additionally, unlike US, FFS uti- 
lizes dynamical simulations and hence this technique does 
not assume that the system is in (quasi-)equilibrium. 

Molecular dynamics and Brownian dynamics (BD) 
simulations are ideal for studying the time evolution of 
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TABLE I: Packing fraction (77 = 7rcr^iV/6V) , reduced pressure 
(/3p(T^) and chemical potential difference between the fluid 
and solid phases( \Ap\) of the state points studied in this pa- 
per. The chemical potential difi^erence was determined using 
thermodynamic integration,'^ and the equations of state for 
the fluid and solid are from Refs. I18I19I respectively. 



systems, and, when possible, they are the natural tech- 
nique to study dynamical processes such as nucleation. 
Unfortunately, however, available computational time of- 
ten limits the types of systems which can be effectively 
studied by these dynamical techniques. Brownian dy- 
namics simulations, which would be the natural choice to 
use for colloidal systems, are very slow due to the small 
time steps required to handle the steep potential used 
to approximate the hard-sphere potential. Event driven 
MD simulations are much more efficient to simulate hard 
spheres and enable us to study spontaneous nucleation 
of hard-sphere mixtures over a range of volume fractions. 
The main difference between the two simulation methods 
regards how they treat the short-time motion of the par- 
ticles. Fortunately, the nucleation rate is only dependent 
on the long-time dynamics which are not sensitive to the 
details of the short-time dynamics of the system.^iSl 

In this paper we study in detail the application of US 
and FFS techniques to crystal nucleation of hard spheres, 
and predict the associated nucleation rates. Combining 
these nucleation rates with results from MD simulations, 
we make predictions for the nucleation rates over a wide 
range of packing fractions ?/ = 0.5214—0.5572, with corre- 
sponding pressures and supersaturations shown in Table 
|lj We compare these theoretical nucleation rates with the 
rates measured experimentally by Refs. HKS 

This paper is organized as follows: in section [Tl] we 
describe and examine the order parameter used to dis- 
tinguish between solid- and fluid-like particles through- 
out this paper, in section [111] we calculate essentially the 
"exact" nucleation rates using MD simulations, in sec- 
tions |IV] and |V] we calculate the nucleation rates of hard 
spheres using US and FFS respectively, and discuss dif- 
ficulties in the application of these techniques, in section 
|VI| we summarize the theoretical results and compare the 
predicted nucleation rates with the measured experimen- 
tal rates of Harland and Van Megen^Sinn et al.^ and 
Schatzel and Ackerson^and section lVlTl contains our con- 
clusions. 
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II. ORDER PARAMETER 

In this paper, an order parameter is used to differenti- 
ate between hquid-hke and soUd-hke particles and a clus- 
ter algorithm is used to identify the solid clusters. For 
this study we have chosen to use the lo cal bo nd-order pa- 
rameter introduced by Ten Wolde et al^^^ in the study 
of crystal nucleation in a Lennard- Jones system. This 
order parameter has been used in many crystal nucle- 
ation studies, including a previous study of hard-sphere 
nucleation by Auer and Frenkel.'* 

In the calculation of the local bond order parameter 
a list of "neighbours" is determined for each particle. 
The neighbours of particle i include all particles within 
a radial distance Tc of particle i, and the total number of 
neighbours is denoted Ni,{i). A bond orientational order 
parameter qi^m{i) for each particle is then defined as 

qi,m{i) = ^ Yl^,n{^^.J,(t>^,3) (3) 

where yj^ml^,'/') are the spherical harmonics, m e [—1,1] 
and Oij and (pi^j are the polar and azimuthal angles of the 
center-of-mass distance vector r.y = Vj — Yi with the 
position vector of particle i. Solid-like particles are iden- 
tified as particles for which the number of connections 
per particle ^(i) is at least and where 

^{i)^Y.H{di{i,3)~d,), (4) 

H is the Heaviside step function, is the dot-product 
cutoff, and 

I 

m 

di[i,J) = Y72 172- 

/ ; \ / ' \ 

\rii— — l / \jn— — l / 

(5) 

A cluster contains all solid-like particles which have a 
solid-like neighbour in the same cluster. Thus each par- 
ticle can be a member of only a single cluster. 

The parameters contained in this algorithm include the 
neighbour cutoff Tc, the dot-product cutoff dc, the criti- 
cal value for the number of solid- like neighbours ^c, and 
the symmetry index for the bond orientational order pa- 
rameter I. The solid nucleus of a hard-sphere crystal is 
expected to have random hexagonal order, thus the sym- 
metry index is chosen to be 6 in all cases in this study. 

To investigate the effect of the choice of ^c, we exam- 
ined the number of correlated bonds per particle at the 
liquid-solid interface. To this end, we constructed a con- 
figuration in the coexistence region in an elongated box 
by attaching a box containing an equilibrated random- 
hexagonal-close-packed (RHCP) crystal to a box contain- 
ing an equilibrated fiuid. Note that the RHCP crystal 



was placed in the box such that the hexagonal layers were 
parallel to the interface. The new box was then equili- 
brated in an NPT MC simulation. We then examined 
the density profile of solid-like particles as determined 
by our order parameter using Tc = 1.4, dc — 0.7 and 

= 5, 7 and 9. As shown in Fig. [T] for all values of 
that we examined the order parameter appears to con- 
sistently identify the particles belonging to the bulk fluid 
and solid regions. For comparison we also show a typical 
configuration of the RHCP crystal in coexistance with 
the fluid phase. The solid-like particles as defined by the 
order parameter are labelled according to the number of 
solid-like neighbours while the fiuid-like particles are de- 
noted by dots. The main difference between these order 
parameters relates to distinguishing between fluid- and 
solid-like particles at the fluid-solid interface. Unsurpris- 
ingly, the location of the interface seems to shift in the 
direction of the bulk solid as is increased. We note that 
the dips in the density profile correspond to HCP stacked 
layers which are more pronounced for higher values of ^c- 
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FIG. 1: Top: A typical configuration of an equilibrated 
random-hexagonal-close-packed (RHCP) crystal in coexis- 
tance with an equilibrated fluid. The crystalline particles are 
labelled according to three different crystallinity criteria: the 
red particle have between ^ = 5 and 6 crystalline bonds, the 
green particles have between ^ = 7 and 8 crystalline bonds 
and the blue particles have ^ > 9 or more crystalline bonds. 
The fluid-like particles < 5) are denoted by dots. Bottom: 
The density profile of particles with a minimum number of 
neighbours ^ as labelled. Note that the dips in the density 
profile correspond to HCP stacked layers. 
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III. MOLECULAR DYNAMICS 

A. Nucleation Rates 

In MD simulations the equations of motion are inte- 
grated to follow the time evolution of the system. Since 
the hard-sphere potential is discontinuous the interac- 
tions only take place when particles collide. Thus the 
particles move in straight lines (ballistic) until they en- 
counter another particle with which they perform an elas- 
tic collision.!^ These collision events are identified and 
handled in order of occurrence using an event driven sim- 
ulation. 

In theory, using an MD simulation to determine nucle- 
ation rates is quite simple. Starting with an equilibrated 
fluid configuration, an MD simulation is used to evolve 
the system until the largest cluster in the system exceeds 
the critical nucleus size. The MD time associated with 
such an event is then measured and averaged over many 
initial configurations. The nucleation rate is given by 



k = 



{t)V 



(6) 



where V is the volume of the system and (t) is the av- 
erage time to form a critical nucleus. Measuring this 
time is relatively easy for low supersaturations where the 
nucleation times are relatively long compared to the nu- 
cleation event itself, which corresponds with a steep in- 
crease in the crystalline fraction of the system. However, 
for high supersaturations pinpointing the time of a nu- 
cleation event is more difficult. Often many nuclei form 
immediately and the critical nucleus sizes must be esti- 
mated from CNT or US simulations. Additionally, the 
precise details of the initial configuration can play a role 
at high supersaturations since the equilibration time of 
the fluid is of the same order of magnitude as the nucle- 
ation time. 

For the results in this paper, we performed MD simu- 
lations with up to 100,000 particles in a cubic box with 
periodic boundary conditions in an NVE ensemble. Time 
was measured in MD units ay^m/ksT. The order pa- 
rameter was measured every 10 time units and when the 
largest cluster exceeded the critical size by 100 percent 
we estimated the time Tnuci at which the critical nucleus 
was formed using stored previous configurations. We per- 
formed up to 20 runs for every density and averaged the 
nucleation times. 

The results are shown in Table HIl The nucleation times 
shown here are for a system of 2.0 • 10* particles and in 
MD time units. To compare with other data we con- 
vert the MD time units to units of /{6Di) with Di the 
long-time diffusion coefficient measured in the same MD 
simulations. We were not able to measure the long-time 
diffusion coefficients for high densities because our mea- 
surements were influenced by crystallization. We used 
the fit obtained by Zaccarelli et who used polydis- 
perse particles to prevent crystallization. For rj < 0.54, 
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TABLE II: The average nucleation time, obtained from MD 
simulations, to form a critical cluster that grew out and filled 
the box. The last column contains the rate (k) in units of 



we find good agreement between our data for and 
this fit. 



IV. UMBRELLA SAMPLING 
A. Gibbs Free-Energy Barriers 

Umbrella sampling is a technique developed by Torrie 
and Valleau to study systems where Boltzmann- weighted 
sampling is inefficient.^^ This method has been applied 
frequently to study rare events, such as nucleation ,1^21 and 
specifically has been applied in the past to study the nu- 
cleation of hard spheres. In general, umbrella sampling 
is used to examine parts of configurational space which 
are unaccessible by traditional schemes, eg. Metropolis 
Monte Carlo simulations. Typically, a biasing potential 
is added to the true interaction potential causing the sys- 
tem to oversample a region of configuration space. The 
biasing potential, however, is added in a manner such 
that is it easy to "un"-bias the measurables. 

In the case of nucleation, while it is simple to sam- 
ple the fluid, crystalline clusters of larger sizes will be 
rare, and as such, impossible to sample on reasonable 
time scales. The typical biasing potential for studying 
nucleation is given bji^nilll 



t/bias(n(r^)) = ^(n(r^)-nc) 



(7) 



where A is a coupling parameter, n{r^) is the size of 
the largest cluster associated with configuration r^, and 
nc is the targeted cluster size. By choosing A carefully, 
the simulation will fiuctuate around the part of config- 
urational space with n(r^) in the vicinity of nc- The 
expectation value of an observable A is then given by 



(A) 



{A/Win{r^))) 



bias 



(l/W^(n(r^))),. 



bias 



where 



W{x) 



(8) 



(9) 



Using this scheme to measure the probability distribution 
P{n) for clusters of size n, the Gibbs free energy barrier 



5 



can be determined bjBSl 

l3AG{n) = constant - ln(P(n)). (10) 

Many more details on this method are given 
elsewhere! ^ "" I ^^ l 




200 



FIG. 2: Gibbs free energy barriers /3AG(n) as a function of 
cluster-size n as obtained from umbrella sampling simulations 
at a reduced pressure of Ppa^ = 17 for varying critical number 
of solid-like neighbours as labelled. For ~ 5, 7, 9, the 
neighbour cutoff is rc — 1.4 and for = 6, 8, 10, = 1.3. In 
all cases the dot product cutoff is dc = 0.7. 
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FIG. 3: Classical nucleation theory fits (thick lines) to the 
Gibbs free energy barriers obtained from umbrella sampling 
simulations at a reduced pressure of /3pa'^ — 17 for varying 
as labelled. Note that the CNT radius (Rcnt) is related to 
the radius (-R(Cc)) measured by umbrella sampling by R{£,c) = 
Rcnt + ^(Cc), where a(^c) is a constant that corrects for 
the different ways the various order parameters identify the 
particles at the fluid-solid interface. The fit parameters are 
given in Table [lV A[ We have shifted the barriers for = 6 — 9 
by 5, 10, 15, 20 fcsT respectively for clarity 
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FIG. 4: Fits of an adjust ed cla ssical nucleation theory 
(ACNT) presented in Section IV A to the Gibbs free energy 
barriers predicted using umbrella sampling simulations at a 
reduced pressure of I3pa^ = 17 and for varying as labelled. 
Note that the CNT radius (i?cNT) is related to the radius 
measured by umbrella sampling by R{^c) = Rcnt + ct{^c), 
where a(^c) is a constant. The fit parameters are given in 
Table |IV A| We have shifted the barriers for = 6 — 9 by 
5, 10, 15, 20 ksT respectively for clarity. 



For a pressure of I3pa^ — 17, corresponding to a super- 
saturation of /3 |A/x| = 0.54, we examine the effect of one 
of the order parameter variables, namely ^c-, on the pre- 
diction of the nucleation barriers. The barriers predicted 
by US using = 5, 6, 7, 8, 9 and 10 are shown in Fig. [2] 
Note that the height of the barriers does not depend on 

within error bars. In general, for larger values of 
more particles are identified as fluid as compared with 
smaller values of ^c- This is consistent with the differ- 
ences between these order parameters as demonstrated 
in Fig. [1] 

Taking the previous discussion on order parameters 
into consideration, we fit the barriers corresponding to 

— 5, 6, 7, 8 and 9 using CNT where we assume there 
exists a CNT radius i?cNT which differs from the radius 
R{^c) measured by the order parameter. We assume that 
the difference (a) is a constant for each value of the crit- 
ical number of solid-like neighbours which corrects for 
the different ways the various order parameters identify 
the particles at the fluid-solid interface: 



i?(Cc) -i?CNT + a(^,). 



(11) 



Note that we have assumed that the cluster size n can 
be related to the cluster radius R{£,c) by 



(12) 



Fitting all barriers simultaneously for the surface ten- 
sion, and the various a{£,c), we obtain the fits displayed 
in Fig. |3] From the various values of a, the associ- 
ated critical CNT radius (i^cNx) 



can be determined. 
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0.54 0.76 


2.49 


-0.425 


-0.231 


-0.000 


0.139 


0.380 








ACNT 


0.54 0.61 


2.01 


-0.961 


-0.765 


-0.551 


-0.402 


-0.148 


8.75 9.46 9.81 


9.78 


9.28 



TABLE III: Numerical values for the parameters associated with the fits in Figs. [3] and [4] for classical nucleation theory and 
the adjusted classical nucleation theory presented in this paper. 



We find R*q^^ = 2.49(t. Additionally, we find a sur- 
face tension of P^a"^ = 0.76 which roughly agrees with 
the results of Auer and Frenkel who obtained surface 
tensions of fi^a'^ = 0.699,0.738 and 0.748 for pressures 
j3pa^ — 15, 16 and 17 respectively.!^ However, recent cal- 
culations by Davidchack et al^^ of the surface tension at 
the fluid-solid coexistence find P^cr"^ = 0.574, 0.557 and 

0. 546 for the crystal planes (100), (110), and (111) re- 
spectively. For a spherical nucleus, the surface tension is 
expected to be an average over the crystal planes. Thus 
our result for the surface tension and that of Ref. 4 ap- 
pear to be an overestimate. 

There have been a number of papers discussing possi- 
ble corrections to CNT (eg. Refs. I27|28p . Recent work 
on the 2d Ising model, a system where both the sur- 
face tension and supersaturation are known analytically, 
demonstrated that in order to match a nucleation barrier 
obtained from US to CNT, two correction terms were re- 
quired, specifically a term proportional to \og{N) as well 
as a constant shift in AG which we define as 0.123 The uS 
barrier is only expected to match CNT near the top of 
the barrier where the log(A^) term is almost a constant. 
Thus, we propose fitting the barrier to an adjusted ex- 
pression for CNT (ACNT), by adding a constant c to Eq. 

1. Fitting the US barriers with this proposed form for the 
Gibbs free energy barrier, where we assume c is a function 
of ^c, we obtain the fits displayed in Fig. |4] In this case 
we find a surface tension P^u"^ = 0.61, and the values for 
a(^c) and c(^c) are given in Table IV A[ The difference in 
the various c(^c) are around Ik^T and correspond well 
to the difference in heights of the barriers. More strik- 
ingly, the surface tension predicted from this proposed 
free energy barrier is in much better agreement with re- 
cent calculations of Davidchack et al. than the surface 
tension we calculate using classical nucleation theory di- 
rectly. We would like to point out here that due to the 
simple form of the nucleation barrier, it is difficult to be 
certain of any fit with more than one fitting parameter, 
as there are many combinations of parameters which fit 
almost equally well. 

Using both expressions for the Gibbs free energy bar- 
rier, namely CNT and ACNT, we were unable to fit the 
barrier corresponding to Ppa^ = 17 and ric — 10 simul- 
taneously with the other predicted barriers for the same 
pressure. We speculate that our difficulty in fitting the 
barrier at = 10 stems from an "over-biasing" of the 
system. Specifically, by using — 10 the biasing po- 
tential could cause the system to sample more frequently 
more ordered clusters, and hence change slightly the re- 
gion of phase space available to the US simulations. In 
general, the least biased systems would be expected to 



explore the largest region of phase space resulting in the 
best results. 

In conclusion, with the exception of — 10, the value 
of used in the order parameter did not appear to have 
an effect on the nucleation barriers once the difference in 
their measurements of the solid-liquid interface was taken 
into consideration. Finally, for use in our nucleation rate 
calculations (section IV B[ ) we also calculated the Gibbs 
free energy AG{n) for reduced pressures f3pa^ — 15 and 
16 using umbrella sampling simulations. We present the 
barrier heights in Table llV) 



B. Umbrella Sampling Nucleation Rates 

The nucleation barriers as obtained from US simula- 
tions can be used to determine the nucleation rates. The 
crystal nucleation rate k is related to the free energy bar- 
rier {G{n)) b>SI 



where 



A^pfn 



\PAG"{n*)\ 



2n 



(13) 



(14) 



n* is the number of particles in the critical nucleus, p 
is the number density of the supersaturated fluid, /„• 
is the rate particles are attached to the critical cluster, 
and G" is the second derivative of the Gibbs free energy 
barrier. Auer and FrenkeP showed that the attachment 
rate /„* could be related to the mean square deviation 
of the cluster size at the top of the barrier by 



1 (An^it)) 



t 



(15) 



The mean square deviation of the cluster size can then 
be calculated by either employing a kinetic MC simula- 
tion or a MD simulation at the top of the barrier. For 
simplicity, in the remainder of this paper the nucleation 
rate determined using this method will be referred to 
as umbrella sampling (US) nucleation rates, although to 
calculate the nucleation rates both US simulations and 
dynamical simulations (KMC or MD) are necessary. 

The mean square deviation, or variance, in the cluster 
size appearing in Eq. [15] has both a short- and long-time 
behaviour. At short times, fluctuations are due to par- 
ticles performing Brownian motion around their average 
positions while the long-time behaviour is caused by rear- 
rangements of particles required for the barrier crossings. 



7 



The slope of the variance is large at short times where 
only the fast rattling is sampled. However, the longer the 
time the further the system has diffused away from the 
critical cluster size at the top of the nucleation barrier. 
AueP^ states that runs need to be selected that remain 
at the top of the barrier. However, when this is done the 
attachment rate is lower than when the average over all 
runs is taken since it excludes the runs that move off the 
barrier fast and have the largest attachment rate. This 
problem is analogous to determining the diffusion con- 
stant of a particle performing a random walk. By only 
including walks which remain in the vicinity of the ori- 
gin, the measurement is biased and excludes trajectories 
which quickly move away from the origin. This results 
is an underestimation of the diffusion constant, and sim- 
ilarly, in this case, an underestimate of the attachment 
rate. In Fig. [5] we demonstrate how, starting from a criti- 
cal cluster, the size of the nucleus fluctuates as a function 
of time and, in fact, can completely disappear or double 
in size within 0.3/t; where r; is the time that it takes a 
particle on average to diffuse over a distance equal to its 
diameter i.e. ti = cr^/(6D;). 
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Time t in IVIC cycles 

FIG. 6: The mean squared deviation (MSD) of the cluster 
size <^An^(f)) as function of time t in MC cycles. The cluster 
size has been measured every cycle and averaged over 100 
cycles to reduce the short-time fluctuations. The slope of this 
graph is twice the attachment rate (Eq. 151. 
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TABLE IV: Nucleation rates k in units of Do/a^ with Do 
the short time diffusion coefficient as a function of reduced 
pressure {Ppa^) as predicted by umbrella sampling. G'^n") 
is the second order derivative of the Gibbs free energy at the 
critical nucleus size n*. 



rates for pressures f3pa'^ — 15, 16, 17 are reported in Table 
lYl 



FIG. 5: The cluster size {n{t)) as a function of time in MC 
cycles for a random selection of clusters that start at the top 
of the nucleation barrier. 

The kinetic prefactor was determined using KMC sim- 
ulations with 3000 particles in an NVT ensemble in a 
cubic box with periodic boundary conditions. The initial 
configurations were taken from US simulations in one of 
the windows at the top of the barrier. We examined 
the results from both Gaussian and normally distributed 
Monte Carlo steps and found agreement within the sta- 
tistical errors. For all the simulations, the MC stepsize 
was between O.Olcr and 0.1a. The variance of the cluster 
size for a typical system is shown in Fig. [6] We observed 
a large variance in the rates calculated for different nu- 
clei. Specifically, some nuclei have attachment rates more 
than an order of magnitude higher than other nuclei of 
similar size. The nuclei with low attachment rates ap- 
peared to have a smoother surface than the nuclei with 
a high attachment rate. 

Our results for the kinetic prefactors and nucleation 



V. FORWARD FLUX SAMPLING 
A. Method 

The forward flux sampling method was introduced by 
Allen et al^^ in 2005 to study rare events and has since 
been applied to a wide variety of systems. Two review 
articles (Refs. 30 31]) on the subject have appeared re- 
cently and provide a thorough overview of the method. 
In the present paper we discuss FFS as it pertains to 
the liquid to solid nucleation process in hard spheres. In 
general, FFS follows the progress of a reaction coordi- 
nate during a rare event. For hard-sphere nucleation, 
a reasonable reaction coordinate (Q) is the number of 
particles in the largest crystalline cluster in the system 
(n). For the remainder of this paper, for all FFS calcu- 
lations, we take the reaction coordinate to be the order 
parameter discussed in Sec. [H] with — Tc ~ 1-3, 
and c?c = 0-7. In general, the reaction coordinate is used 
to divide phase space by a sequence of interfaces (Aq, 
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Ai, ... Ajv) associated with increasing values n(r^) such 
that the nucleation process between any two interfaces 
can be examined. In our case the hquid is composed of 
all states with n < Ao and the solid contains all states 
with n > Aat. While the complete nucleation event is 
rare, the interfaces are chosen such that the part of the 
nucleation process between consecutive interfaces is not 
rare, and can thus be thoroughly studied. 

In the FFS methodology, the nucleation rate from the 
fluid phase A to the solid phase B is given by 



kAB = $AAo^'(AAr|Ao) 



(16) 
(17) 



1=0 



where ^aXq is the steady-state flux of trajectories leaving 
the A state and crossing the interface Xq in a volume V, 
and P(Ai+i|Ai) is the probability that a configuration 
starting at interface A^ will reach interface A^+i before it 
returns to the fluid {A). 




500 1 000 

Time t in IVIC cycles 



1500 



FIG. 7: The cluster size as a function of time t in MC cycles 
for 4 random trajectories at pressure Ppa"^ — 17 starting with 
a cluster size of n = 43 using kinetic MC simulations with 
stepsize Akmc ~ 0.1a and measuring the order parameter 
every Atord = 5 MC steps. 

If we apply this method directly to a hard-sphere sys- 
tem a number of difhculties arise. As shown in Fig. [5j on 
short times the size of a cluster measured by the order 
parameter fluctuates wildly. The variance in the cluster 
size displays two different types of behaviour, short-time 
fluctuations related to surface fluctuations of the clus- 
ter, and a longer time cluster growth (Fig. [6|. Thus, if 
we try to measure the flux $aAo directly, we encounter 
difficulties due to these short-time surface fluctuations. 
In theory, FFS should be able to handle these types of 
fluctuations, however, they increase the amount of statis- 
tics necessary to properly measure the flux and the first 
probability window properly. In the second part of FFS 
calculations, probabilities of the form P(Ai+i|Ai) need to 



be determined. In calculating these probabilities it is im- 
portant to be able to determine if a cluster has returned 
to the fluid (A). For pre-critical clusters we find large 
fiuctuations of the order parameter, as shown in Fig. [Tj 
which can lead to a cluster being misidentified as the 
fiuid (A). Specifically, in this figure the darkest trajec- 
tory (black) shows a cluster containing 43 particles that 
shrinks to 5 particles before it returns to 40, and finally 
reaches a cluster size of 60 particles. Hence, if we had 
set Ao = 5, this trajectory would have been identified as 
melting back to the fluid phase (A). However, since the 
growth of a cluster from size 5 to 60 is a rare event in 
our system, we presume that this was simply a short- 
time fluctuation of the cluster and not a 'real' melting 
of the instantaneously measured cluster. For pre-critical 
clusters, these fluctuations result in cluster sizes that are 
smaller than the cluster 'really' is. We suggest that these 
fluctuations are largely related to the difficulty that this 
order parameter has in distinguishing between solid- and 
fiuid-like particles at the fluid-solid interface. For larger 
clusters, where the surface to volume ratio is small, this 
problem is minimal. However, for elongated or rough 
pre-critical clusters, where the surface to volume ratio 
is large, these surface fluctuations and rearrangements 
are important, and can cause problems in measuring the 
order parameter. 

Thus, to try and address these problems, in this paper, 
we apply forward flux sampling in a slightly novel way. 
We regroup the elements of the rate calculation such that 



kAB = $AAi n 



i=l 



where 



«'AAi = $AAo^(Al|Ao). 



(18) 



(19) 



We note that if Ai is chosen such it is a relatively rare 
event for trajectories starting in A to reach Ai, then 



1 



A\i 



(20) 



where (tAXi) is the average time it takes a trajectory in A 
to reach Ai. The approximation made here, in contrast 
to normal FFS simulations, is that the time the system 
spends with an order parameter greater than Ai is neg- 
ligible. Since even reaching this interface is a rare event, 
this approximation should have a minimal effect on the 
resulting rate. Additionally, in this way we are relatively 
free to place the flrst interface (Ao) anywhere under XiW^ 
We choose to use Aq = 1 to minimize the effect of fluctu- 
ations, as seen in Fig. [Tj on the probability to reach the 
following interface. Here we assume that any crystalline 
order in a system with an order parameter of 1 likely does 
not arise from fluctuation of a much larger cluster, but 
rather is very close to the fluid, and is expected to fully 
melt and not grow out to the next interface. In this man- 
ner we are able to start several parallel trajectories from 
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the fluid in order to measure {Iaxi), stopping whenever 
the trajectory first hits interface Ai. 

In our implementation of FFS, we employ kinetic 
Monte Carlo (KMC) simulations at fixed pressure to fol- 
low the trajectories from the liquid to the solid. The 
KMC simulations are characterized by two parameters, 
the maximum stepsize (Akmc) per attempt to move each 
particle, and the frequency with which the order param- 
eter (reaction coordinate) is measured A^ord- However, 
during an FFS simulation, it is expected that the order 
parameter is known at all times such that it is possible to 
identify exactly when and if a given simulation reaches 
an interface. Thus it is possible that A^ord introduces an 
additional error into our measurement of the rate. 

To examine the effects of i) the approximation associ- 
ated with our method for calculating I'aAi , ii) the short- 
time fluctuations of the order parameter (which could be 
considered as an error in the measurement of the clus- 
ter size), and iii) the frequency of measuring the order 
parameter, we examined the nucleation rate for a simple 
one-dimensional model system in the presence of such 
features. Details of these simulations are given in Ap- 
pendix A. In this simple model system, we find that none 
of these features have a large effect on the rate. In fact, 
for most cases, the difference is too small to see within 
our error bars. 



B. Simulation details and results 

All simulations were performed with 3000 particle in a 
cubic box with periodic boundary conditions. Initial con- 
figurations were produced using NPT MC simulations of 
a liquid phase at a reduced pressure of /3pa^ — 1000. The 
simulations were stopped when the packing fraction as- 
sociated with the pressure of interest was reached. This 
initial configuration was then relaxed using an NPT sim- 
ulation at the correct pressure {f3pa^ = 15, 16, 17). The 
relaxation consisted of at least 10,000 MC cycles, after 
which the simulation continued until a measurement of 
the order parameter found no crystalline particles in the 
system. 

In order to determine the flux and the probabilities, 
100 trajectories were started in the liquid and termi- 
nated when n{r-^) — Ai. These trajectories were pro- 
duced using KMC simulations. The probability P(A2|Ai) 
was then found by making C'l copies of the configurations 
that reached Ai, and following these configurations until 
they either reached A2 or returned to the fluid. By tak- 
ing different random number seeds, the various copies of 
the same conflgurations follow different trajectories. The 
fraction of successful trajectories corresponds to the re- 
quired probability. The successful trajectories were then 
copied C2 times to determine P(A3|A2). The remaining 
P{K+i\^i)'s are calculated similarly. 

To study the effect of the two KMC parameters, 
namely Akmc £^nd Aiordj on the nucleation rates, we 
have examined the first 8 FFS windows for l3pa^ = 15 



for various values of the number of MC steps between the 
order parameter measurements Aiord and the maximum 
displacement Akmc for the KMC simulations. The re- 
sults are shown in Table |Vl As shown in this table we do 
not find a significant effect on the rate from either pa- 
rameter. Thus for numerical efficiency, unless otherwise 
indicated, the rates in this section come from Aiord = 5 
MC cycles and Akmc = 0.2a. 

For pressures Ppa^ = 16 and 17 we have performed two 
separate FFS calculations to determine the nucleation 
rates, and for pressure I3pa^ = 15 we have the result 
from a single FFS simulation. A summary of the results 



are given in Table VI A complete summary of the results 
for P(Ai-|-l|Ai) for each simulation is given in Tables VII 
[Villi and[lXl 



VI. SUMMARY AND DISCUSSION 
A. Nucleation Rates 

In this section we examine hard-sphere nucleation rates 
predicted using US simulations, MD simulations and 
FFS simulations together with the experimental results 
of Harland and Van Megen,!!] Sinn et al.^ and Schatzel 
and AckersorP and the US simulations of monodisperse 
and 5% polydisperse hard-spheres mixtures examined 
by Auer and Frenkel.^I The experimental volume frac- 
tions have been scaled to yield the coexistence densities 
of monodisperse hard spheres.'i^ Similarly, we scale the 
polydisperse results of Auer and Frenkel with the coex- 
istence densities determined in Ref. [3H Inspired by the 
recent work of Pusey et al.^we plot the nucleation rates 
in units of the long-time diffusion coefficient. In experi- 
ments with colloidal particles, the influence of the solvent 
on the dynamics cannot be ignored. Specifically, the sys- 
tem slows down due to hydrodynamic interactions when 
the density is increased. However, since hydrodynamics 
are included in the long-time diffusion units, if we present 
the nucleation rates in terms of the long-time diffusion 
coefficient, our predicted nucleation rates should be in 
agreement with the experiments. The time in experi- 
ments is typically measured in units of Dq, the free diffu- 
sion at low density. We convert the short-time diffusion 
coefficient Dq to long-time diffusion coefficient Di using 



Do 



V 0.58/ 



(21) 



Harland and Van MegerF* claim that S ~ 2.6 gives a good 
fit to their system and Sinn et al.^ use S — 2.58. Since the 
system Schatzel and AckersorP examine is very similar to 
the other two, we use S = 2.6 to convert their nucleation 
rates to long-time units. We note that both S — 2.58 and 
6 = 2.6 give very similar results. The results for both the 
theoretical and experimental rates in long time units are 
shown in Fig. [H] 

In Ref. [T^ Pusey et ai showed that the nucleation 
rates for various polydispersies (0 to 6%) of hard-sphere 
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TABLE V: Probabilities P{\i+i\Xi) for the first 8 interfaces for a pressure of /3pa^ = 15 where the KMC simulations stepsize 
(Akmc) and the number of MC steps between measuring the order parameter Atord are varied. The following interfaces were 
used: A2 = 20, A3 — 26, A4 = 32, A5 = 38, Ag — 44, A7 = 54, As = 65, and Ag = 78. In all cases, 100 configurations were started 
in the fluid and reached the first interface, and at each interface, d — 10 copies of each successful configuration were used. 
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FIG. 8: A comparison of the crystal nucleation rates of hard spheres as determined by the three methods described in this 
paper FFS, US, and MD with the experimental results from Refs. [JHS] and previous theoretical results from Ref. 3j Note 
that error bars have not been included in this plot. In general, the error bars of the simulated nucleation rates are largest 
for lower supersaturations (ie. lower volume fractions), as the barrier height is higher. For the FFS and US simulations, the 
error for Ppa^ = 15 (77 = 0.5214) is between 2 and 3 orders of magnitude, and for /3pcr^ = 17 (r; = 0.5352) is approximately 
one to two orders of magnitude. The MD results are quite accurate around /3pcr^ = 17, however the error bars are larger for 
the higher pressure MD results. Within these estimated error bars, the simulated nucleation rates are all in agreement, while 
the experimentally obtained rates show a markedly different behaviour, particularly for low supersaturations where the the 
difference between the simulations and experiments can be as large as 12 orders of magnitude. 



mixtures collapsed onto the same curve when the rates 
were plotted in units of the long-time diffusion coefScient. 
We find similar results here. Both the monodisperse and 
polydisperse US results of Auer and Frenkel,* in addition 
to our own US predictions of the nucleation rate, agree 
well within the expected measurement error. Addition- 
ally, we find that the simulation results of the US, FFS, 
and MD all agree. 

However, on the experimental side, the nucleation rates 



of Harland and Van MegerPare approximately one to two 
orders of magnitude below the experiments of Sinn et al^ 
and Schatzel and Ackerson.l^ This is unexpected due to 
the similarity between the experimental systems. In our 
opinion, the main difference between these experiments 
is the polydispersity of the particle mixtures: 5% in the 
case of Harland and Van Megen,'i' 2.5% in the case of 
Sinn et al.,^ and < 5% for Schatzel and Ackerson.'^ How- 
ever, as demonstrated by Pusey et aZ.,— and now also in 
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TABLE VL Nucleation rates predicted using forward flux 
sampling in short-time diffusion coefficient units (-Do)- The 
probabihties P(A_b!Ai), number of steps between the order 
parameter measurements Aord, and kinetic MC stepsize are 
as in Tables |VII[ |VIII[ and |IX[ At each interface, d copied 
of each successful configuration were used. 
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TABLE VIL Probabilities P(Ai+i|Ai) for the interfaces used 
in calculating the nucleation rate for pressure Ppa^ = 17 with 
step size Akmc = O.lcr and measuring the order parameter 
every Atord = 5 MC cycles. 
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Fig. [Sj the nucleation rate when measured in long-time 
diffusion coefficient units should not be effected by the 
polydispersity. Thus, this seems unlikely as an explana- 
tion. A more probable difference between the results may 
simply be due to measurement error in the experimental 
volume fraction which is extremely difficult to measure. 
As shown in Fig. [8] a measurement error of Ary = 0.05 
can have a very large effect on the nucleation rates. It 
is worth mentioning here that the units of Ref . [T] were 
not always mentioned, and thus there also may be some 
error in the manner in which we converted these rates to 
long-time units. 

When we compare the experimental rates with the the- 
oretical results, we find that while the experiments ap- 
pear to match the general trend of the simulations for 
high supersaturations they predict a significantly higher 
nucleation rate at lower densities. The argument pre- 
sented above regarding measurement error in the volume 
fraction is one possible explanation, ie. by simply shifting 
the experimentally predicted nucleation rates to higher 
densities the agreement is much better. However, we 
speculate that there is another possible reason for the 
discrepancy. Specifically, at high supersaturations there 
should be many nucleation events occurring in the ex- 
perimental system during a fairly short time interval. 
Hence, it should be possible to measure the nucleation 
rate before a single cluster has the chance to grow out 
significantly. However, at lower supersaturations, when 
a nucleation event is extremely rare, a single cluster in 
the experimental system can grow out significantly before 
sufficient nucleation events have occurred to measure the 
rate. However, these large clusters can contain a number 
of twinning defects,!^ resulting in scattering from the var- 
ious crystalline domains. Scattering from these domains 
may lead to an over-count in the nucleation events per 
volume and time unit, yielding higher nucleation rates. 



TABLE VIIL Same as Table VII but for Ppa 
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TABLE IX: Same as Table |VII| but for /3pcr'^ = 15 and with 

Atord ~ 2. 



To examine whether the structure and shape of the 
critical clusters from US simulations depended on the 
precise threshold values used for the crystalline order pa- 
rameters, we compared and analysed the critical clus- 
ters obtained when three different crystalline order pa- 
rameters were used to bias the US simulations, namely, 

= 5, 7 and 9. Subsequently we analyzed these critical 
clusters using the three different order parameters. In 
Fig. [9j two typical critical clusters from different bias- 
ing order parameters are shown on the top and bottom 
rows. The nucleus of the cluster, shown in blue, was 
identified by all three cluster criteria (^c = 5,7 and 9). 
The main difference between the criteria is the location 
of the fluid-solid interface as shown by the green and red 
particles. The strictest order parameter finds only the 
more ordered center whereas the loosest version detects 
the more disordered particles at the interface as well. 



If Fig. 10 we show some of the nuclei obtained from 
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ec = 5 = 7 ^'^T^ 




FIG. 9: Two typical snapshots (top and bottom) of the critical nuclei as obtained with US at a volume fraction rj — 0.5355 
using different values of the critical number of crystalline bonds = 5 (left), 7 (middle) and 9 (right) in the biasing potential. 
The clusters are analyzed with three different crystalline order parameters. The blue particles are found by all three cluster 
criteria, the green particles have ^ = 7 or 8 crystalline bonds and the red particles have only ^ = 5 or 6 crystalline bonds. 



MD simulations. These snapshots were taken just before 
the nuclei grew out so they are not necessarily precisely 
at the top of the nucleation barrier. They appear very 
similar in roughness and aspect ratio to those obtained 
from US simulations. 

To further examine whether the choice of method in- 
fluenced the resulting clusters, we calculated the radius 
of gyration tensor for each of the methods for pres sur e 
Ppa'^ = 17 as a function of cluster size (see Figure 11). 
There is no indication that the clusters in any of the sim- 
ulation methods differed substantially. 

Additionally, we examined whether the simulation 
technique influenced the type of pre-critical nuclei that 
formed in the simulations, ie. face-centered-cubic (FCC), 
and hexagonal-close-packed (HCP). To do this we used 
the order parameter introduced by Ref. 34 which allows 
us to identify each particle in the cluster as either FCC- 
like or HCP-like. The results for a wide range in nucleus 
size is shown in Fig 
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We find complete agreement 
between the three simulation techniques. Specifically, in 
all cases we find that the nucleus is composed of approx- 
imately 80% FCC-like particles. This was unexpected 



as the free energy difference between the bulk FCC and 
HCP phases is about O.OOlksT per particle at meltinJ^H 
and hence a random stacking of hexagonal layers in the 
nuclei would be expected.l^Sl We speculate that this pre- 
dominance of FCC stacking in the nuclei arises from sur- 
face effects. 



VII. CONCLUSIONS 

In conclusion, we have examined crystal nucleation of 
hard spheres with molecular dynamics, umbrella sam- 
pling and forward flux sampling simulations. We find 
that the nucleation rates predicted by all three methods 
agree over the large range in volume fractions we exam- 
ined. Additionally, in agreement with the recent work 
of Pusey et al. we find that by measuring the nucle- 
ation rates in terms of the long-time diffusion constant 
and scaling to the coexistence density of monodisperse 
hard spheres, the 5% polydisperse results of Auer and 
Frenkel'* also agree. On examining the critical clusters, 
we do not find a difference in the nuclei formed using the 
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FIG. 10: Snapshots of spontaneously formed nuclei during an MD simulation at a volume fraction of r\ - 
were taken just before the nuclei grew. The color coding of the particles is the same as in Fig. [9] 
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FIG. 11: A comparison of the three components of the radius 
of gyration tensor as a function of cluster size n, as well as 
the sum of the three components, for clusters produced using 
FFS, MD, and US simulations. 



three simulation techniques. 

We have also compared our nucleation rates with pre- 
vious experimental data, specifically, the nucleation rates 
predicted by Harland and Van MegenjE Sinn et ai^ and 
Schatzel and Ackerson.l^l The nucleation rates measured 
by these three experiments, in contrast to what would 
be expected, differ by about one order of magnitude. 
In general, the experimental systems are similar enough 
that one would have expected agreement in the rate once 
the rates were scaled to the coexistence densities of hard 
spheres. Additionally, while the simulation results agree 
well with the experimental results for high supersatura- 
tions, there is a significant difference between the simula- 
tions and experiments for smaller volume fractions. We 
speculate here that this difference may be due to difficul- 
ties in distinguishing between separate nuclei domains in 
the experiments, or measurement error in the experimen- 



FIG. 12: Fraction of particles identified as either FCC or 
HCP respectively in the clusters produced via molecular dy- 
namics (MD), forward flux sampling (FFS), and umbrella 
sampling (US) simulations as a function of cluster size n. All 
three methods agree and find the pre-critial clusters prodom- 
inately FCC. 



tal volume fractions. 
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TABLE X: Nucleation rates for the one-dimensional potential given by Eq. |A1| and shown in Fig. [13] for Atoid as indicated. 
For each Atord , we performed 10 independent EES simulations. The average rate and associated standard deviation is also 
as indicated. In all cases, 100 configurations were started in the fluid, and at each interface d = 10 copies of the successful 
configurations were used to calculate the proceeding probabilities. The interfaces were placed at Ao = 0, Ai = 1.5, A2 = 1.7, 
A3 = 1.9, A4 — 2.2, A5 = 2.6, Ag — 3.3, and A7 — 4.0 and the flux was calculated using Eq. |20| 
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TABLE XI: Nucleation rates for the one-dimensional potential given by Eq. |Al| and shown in Eig. [Tsjwhere the order parameter 
is given by Eq. |A2| and (JCauss is as indicated. For each (JGauss, we performed 10 independent FES simulations. The average 
rate and associated standard deviation is also as indicated. In all cases, 100 configurations were started in the fluid, and at 
each interface d — 10 copies of the successful configurations were used to calculate the proceeding probabilities. The interfaces 
were placed at Ao = 0, Ai = 1.5, A2 — 1.7, A3 = 1.9, A4 = 2.2, A5 — 2.6, Ae — 3.3, and A7 — 4.0 and the flux was calculated 
using Eq. |20] 



Appendix A: FFS in tiie presence of measurement 
error 

As mentioned in Section|V]of this paper, the FFS tech- 
nique assumes that the reaction coordinate is known ex- 
actly at all times. However, for the hard-sphere system 
examined in this paper, this is not possible due to the 
computational time required for measuring the order pa- 
rameter. In applying the FFS technique to hard spheres, 
two separate types of error are introduced: i) error asso- 
ciated with our inability to know the value of the reaction 
coordinate at all times, and ii) an error in measuring the 
number of particles in a cluster for a given configuration. 
Additionally, as discussed in Section |Vj in this paper we 
have applied FFS in a slightly novel manner. In this ap- 
pendix, we introduce a simple model to examine the ef- 
fect this approximation and the effect such measurement 
errors have on the nucleation rate predicted by forward 



flux sampling. 

To this end, we study the transition rate for a single 
Brownian particle to surmount a one dimensional poten- 
tial energy barrier given by 

l3U{x) ^8x^ -2x^. (Al) 

A plot of the barrier is shown in Fig. [13) For this poten- 
tial, we consider the 'liquid' state to be near x = and 
the 'solid' phase to be near x = A. 

We first determine the 'exact' nucleation rate using 
spontaneous simulations. To do this we perform a ran- 
dom walk starting aX x = and determine the time it 
takes the random walk to surmount the barrier. The 
rate is then given by i? = 1/ (^) • Performing 40 such ran- 
dom walks we find the nucleation rate to be 1.5 • 10~^^. 
In all the calculations in this section, we set the KMC 
stepsize equal to Akmc = 0.1. 

Secondly we explore the effect on the nucleation rate 
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FIG. 13: Toy model potential used to study forward flux 
sampling in the present of various types of measurement error. 



of not knowing the value of the order parameter at all 
times. For this purpose we have performed FFS simu- 
lations when the order parameter was measured every 
Aioid — 1, 2, 5, 10, 50 kinetic Monte Carlo steps. The re- 
sults are shown in Table |X] The average nucleation rates 
predicted for all values of Aiord clearly are the same 
within error. Similarly, the standard error associated 
with Atord = 1, 2, 5, 10 are approximately the same, and 
is only marginally larger for Aiord = 50. Hence we con- 
clude that the frequency of measuring the order parame- 



ter does not significantly affect the predicted nucleation 
rate. Additionally, these nucleation rates agree with the 
nucleation rate predicted from spontaneous simulations 
indicating that of applying FFS as outlined in Section [V| 
predicts the correct nucleation rates. 

Finally, we examine the effect that the measurement 
error in the cluster size has on the nucleation rate. For 
this purpose, we apply a noise term to our order param- 
eter such that 



(A2) 



where Xm is the value of the order parameter used in the 
FFS simulation, cctrue is the true value of the order pa- 
rameter, and S is taken from a Gaussian distribution with 
a mean of and a standard deviation (TGauss- In Table 
|XI| we demonstrate the effect on the predicted nucleation 
rate for various choices of acauss- The resulting nucle- 
ation rates are in good agreement with the spontaneous 
results. For larger CGauss, eg. tJoauss = 0.08 and 0.1, the 
standard error in the results is slightly larger, however, 
the predicted nucleation rates are still correct. 

In summary, we have examined the effect of the ap- 
proximation described by Eq. 20 as well as the effect of 



measurement error in the order parameter and the mea- 
surement frequency Atord of the order parameter. We do 
not find a significant effect on the predicted nucleation 
rates. Thus we conclude that FFS should be robust to 
the types of error we are introducing when we apply the 
technique to hard spheres. 
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